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Abstract 



In nonparametric statistical problems, we wish to find an estima- 
tor of an unknown function /. We can split its error into bias and 
variance terms; Smirnov, Bickel and Rosenblatt have shown that, for 
a histogram or kernel estimate, the supremum norm of the variance 
term is asymptotically distributed as a Gumbel random variable. In 
the following, we prove a version of this result for estimators using 
compactly-supported wavelets, a popular tool in nonparametric statis- 
tics. Our result relies on an assumption on the nature of the wavelet, 
which must be verified by provably-good numerical approximations. 
We verify our assumption for Daubechies wavelets and symlets, with 
A'^ = 6, . . . , 20 vanishing moments; larger values of N, and other wavelet 
bases, are easily checked, and we conjecture that our assumption holds 
also in those cases. 



1 Introduction 

In nonparametric statistical problems, such as density estimation, regres- 
sion, or white noise, we wish to find an estimate / of an unknown function 
/ (Tsybakov, 2009). We can measure the accuracy of an estimator / by its 
distance from /, ||/ — /||, where || • || is some norm on functions. We can 
then decompose the error into variance and bias terms, 

11/ -/II < II/-E/II + IIE/-/II, 

where the bias term ||E/— /|| is deterministic, and the variance term ||/— E/|| 
we hope has an asymptotic distribution independent of /. 
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In density estimation, for the supremum norm on [0, 1], 

||/IL:= sup |/(x)|, 

a:G[0,l] 

the hmiting distribution of a suitably scaled variance term is given by 
Smirnov (1950) for histograms, and in the classical paper of Bickel and 
Rosenblatt (1973) for kernel estimates. In both cases, as the sample size n 
tends to infinity, the variance term approaches a Gumbel distribution, 

oo / / 

for known sequences An, Bn- This result has been of key importance for a 
variety of problems in nonparametric statistics. 

Wavelets are an increasingly popular statistical tool, allowing a simple 
theoretical description of nonparametric problems, and a computationally 
efficient implementation of their solution. Gine and Nickl (2010) establish an 
equivalent of these Smirnov-Bickel-Rosenblatt theorems for certain wavelet 
estimators, using a result of Hiisler, Piterbarg, and Seleznjev (2003) on the 
convergence of cyclostationary Gaussian processes. Gine and Nickl describe 
the asymptotic distribution of the supremum, on increasing intervals, of the 
Gaussian process 

X{x) := j K{x,t)dBt, 

where K is a wavelet projection kernel, and B a Brownian motion; they 
then link this result to the statistical problem considered above. 

Their result holds only for wavelets satisfying certain analytic conditions, 
which the authors demonstrate are satisfied by Battle-Lemarie wavelets hav- 
ing < 4 vanishing moments; Gine, Giintiirk, and Madych (2011) extend 
this to larger values of N. Past work has not, however, succeeded in estab- 
lishing results for the most commonly used wavelets, such as Daubechies 
wavelets and symlets. These wavelets, unlike those of Battle and Lemarie, 
are compactly supported, allowing the most efficient implementation of sta- 
tistical procedures. In the following, we demonstrate that the conditions of 
Gine and Nickl (2010) hold also in these cases, thereby proving a Smirnov- 
Bickel-Rosenblatt theorem for the most practically relevant wavelet bases. 

We work primarily in the white noise model, but also discuss conse- 
quences for the density estimation and regression models. We consider 
wavelet bases both on M, and also on the interval, using the construction of 
Cohen et al. (1993). In both cases, we show that the variance term again 
approaches a Gumbel distribution. We also extend a theorem of Hiisler et al. 
(2003) (as reported in Hiisler, 1999), establishing a uniform convergene re- 
sult for cyclostationary processes; this allows us to show that convergence 
to Gumbel occurs uniformly in large values of the level x. These results are 
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used in Bull (2011) to construct adaptive confidence bands for nonparamet- 
ric statistical problems, and are also of relevance to many other wavelet 
procedures. 

To prove our results, we must first verify an assumption on the wavelet 
functions, which in general do not have an analytic form. We therefore make 
use of provably-accurate numerical approximations, given by Rioul (1992); 
these approximations also provide an efficient means of computing the con- 
stants in our results. We verify our assumption for Daubechies wavelets and 
symlets, having = 6, . . . , 20 vanishing moments (Daubechies, 1992, §6.4); 
however, the numerical approximations can easily be applied to larger values 
of N, and other wavelet bases, and we conjecture that our assumption holds 
also in these cases. 

We state our result in Section 2, and describe the necessary numerical 
approximations in Section 3. We give proofs in Appendix A, and source 
code in Appendix B. 

2 Results 

To begin, we will need (p and ^|J, the scaling function and wavelet of an 
orthonormal multiresolution analysis on L^(M). (For an introduction to 
wavelets and their statistical applications, see Hardle et al., 1998.) We 
make the following assumptions on ip and ip, which are satisfied, for exam- 
ple, by Daubechies wavelets and symlets, with > 6 vanishing moments 
(Daubechies, 1992, §6.1; Rioul, 1992, §14). 

Assumption 2.1. 

(i) For K gN, if and ^ are supported on the interval [1 — K,K]. 
(a) For N gN, ip has N vanishing moments: 



(Hi) if is twice continuously differentiable. 

We will consider wavelet bases on both R and [0, 1], constructed from if 
and ijj. On M, we have an orthonormal basis of L^(R) given by 



for some lower resolution level jo S Z, j > jo, and G Z. 

On [0,1], we can generate an orthonormal basis of L^([0, 1]) using the 
construction of Cohen et al. (1993) (see also Chyzak et al., 2001). We obtain 
basis functions 





ipj,,k{x) := 2^«/V(2^°x - A:) 



iPj^k ■■= 2^'/V(2-'x - k) 



(2.1) 



k = 



...,2^0-1, 
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and 

^j,k, j > Jo, = 0, . . . ,2^ - 1. 

For k € [N, 2^ —N), these functions are given by (2.1). For other values of k, 
the basis functions are specially constructed, so as to form an orthonormal 
basis with desired smoothness properties. 

We will also need to make an assumption on the precise form of the 
scaling function (p. While this assumption is difficult to verify analytically, 
we will see in the following section it can be tested using provably good 
numerical approximations. 

Assumption 2.2. The 1-periodic function 

attains its maximum at a unique point S [0, 1), and (o'^)"(to) < 0. 

Given these assumptions, suppose we have an unknown function /, with 
empirical wavelet coefficients afc,/3j,fc, 

k j>jo k 

Suppose also that we observe the empirical wavelet coefficients 

Ofc := Ofc + ejo,fc' (^j,k '■= l^j,k + ^j,k, (2.2) 

where the t are i.i.d. N{0,a^). This is the case in the white noise model, 
where we observe the process 

Yt= f f{s)ds + n-^'^Bt, 
Jo 

for a Brownian motion B. The empirical wavelet coefficients 
Ofc = y fjoA*) dYt, Pj^k = j ipj,k{*) dYt, 

satisfy (2.2) with = n'^ (Hiirdle et al., 1998, §10). The model (2.2) 
also serves as a limiting approximation in density estimation and regression, 
which we return to later. 

The wavelet projection estimate of /, at resolution level j, is then 

k io<'<i k 

Set 

^ ._ Y:k^i.^'{to-k? 
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and define the quantities 



aij) := V21og(2)j 
b{j) := a{j) 



log(7rlog2) +logi - ilog(l + 



2a(j) 
c{j) := ^2^V2^ 
x(7) := -log(-log(l -7)) . 
We then have the following result on the distribution of the variance term. 
Theorem 2.3. Let jn — s- 00, 70 G (0, 1), and either: 
(i) for a wavelet basis on M, r„ := (0,70]; or 

(a) for a wavelet basis on [0,1], r„ := [7^,70]) where 7„ S (0,70), and 
Tn"*^ = o(e'^-'") for any C > 0. 

Then, as n —)• 00, 

/ x(^) \ \ 

^ 0. 



sup 

7Gr„ 



7-'F (ll/(in) - mJn)\L > C{jn) + 



While this result is stated for the white noise model, similar results hold 
also in density estimation and regression. In density estimation, / is a 
density, and we observe 

Ai, . . . , A„ ~ /. 

This can be linked to the white noise model using Gine and Nickl (2010, 
§4.1). In regression, we have independent observations 

Yi^N{f{xi),a'), 

for Xi := i/n, i = 1, . . . ,n. Regression is known to be asymptotically equiv- 
alent to white noise, as in Brown and Low (1996). We can thus transfer our 
result also to these models. 



3 Numerical approximations 

To apply our result, we must first verify Assumption 2.2, which depends on 
the function if and its derivatives. In general, 99 has no explicit form, but 
we can approximate it numerically using the cascade algorithm, ip satisfies 
a two-scale relation, 

2K-1 

^{x) = uP^{2x + K -k), 

k=0 
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for filter coefficients u'^\ . . . , ^ ^ satisfying 



E 4°' = E 4°' = 1. 

k odd fc even 



and we can use these filter coefficients to compute an approximation to (/?. 
Forn > 0, A; = 0, . . . , 2K - 1, set 



i=0 

and for j,n>0,0<k< T{2K - i), x G [1 - K, K), 
An) ._ /^^^V-lVo^"^ f^"V2;) — 

i=0 ^ ^ 



The functions /j*^^ then converge to a limit function / defined by the ^ 

and the /j'^'' likewise converge to /^"^. The following theorem bounds the 
error in this approximation, and is a straightforward consequence of results 
in Rioul (1992). 

Theorem 3.1. For integers j,n > 0, set 

/ 2K~2 \ 



(n) , 2^-1 I (n+1) 

' — 1 - J log2 I max 2_^^"' 

\ ~ i=0 



a. -.= 1-3 iog2 max 2^ I ' 



Y 1=0 k=0 



max 

m=0, 
V fc=0 



k 

,(n) 



i=0 



If af > for some j, the functions /j"^ converge in L°° to a function 
/ : [1 - K) ^ M satisfying 



2K-1 

,(0) 



f{x)= uff{2x + K-k). 



k=0 

If also aj^^ > for some j, and n > 0, then f is n -times- differentiable, 

and the /j"^ converge in L°° to /("). For n > 0, the approximations /j"^ 
converge at a rate 



(n) 
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Furthermore, given integers j > 0, a < b, set 

I :=2-^[a,b+l) + Z, J{1) ■= \[2^'^ a\ - 2K + 2, [2^'^b\ \ +2^Z. 

Then, on In [1 - K,K): 

(i) for I > j, the values of fj^""^ depend on g^^^ only for k G and 

(a) the above results hold also for quantities a^j^\l) and Cj'^\l) defined 
similarly, taking maxima over g^^^^^^ and fi^^^^ only for k £ J{1)- 

The function ip may be defined as tlie limit of this procedure (Daubechies, 
1992, §6.5). We may thus compute upper and lower bounds on ip and its 
derivatives. Note that, while we could obtain values of the derivatives by 
finite differencing, this would be numerically unstable, and lead to poor 
bounds; the above procedure provides good bounds on all derivatives of ^p. 

To verify Assumption 2.2, and to compute the constants and u<^, we 
must use these bounds to control the function cr^, and its derivatives. Doing 
so over the whole of [0, 1] requires memory exponential in j, which quickly 
becomes infeasible. However, once we have approximated cr^ well enough to 
know that its maxima lie in some interval I, we can exploit the local nature 
of the cascade algorithm, and its bounds, to approximate o"^ only over /. 
As the resolution j increases, so does the accuracy with which we can locate 
the maxima, ensuring our memory costs remain manageable. 

In our implementation, we choose / to be the smallest interval containing 
all points t for which the bounds on cj^, and its derivative, are consistent 
with: 

(i) '^lit) = suPse[o,i] f^^(s); and 



(ii) = 0. 



Note that to ensure efficiency, we must allow choices of / which wrap around 
the edges of [0, 1]; in other words, we must allow I to be any interval on the 
torus. If we find an interval / containing all maxima of a^, with the property 
that C7^ < — e < on /, we may conclude Assumption 2.2 is satisfied. We 
have thus described Algorithm 1. 

To obtain high accuracy, the computation of the filter coefficients Ufc, 
and subsequent approximations, must be performed using variable-precision 
arithmetic; the rounding error in these computations must likewise be con- 
trolled with interval arithmetic. We satisfy these requirements by imple- 
menting the above algorithm in the computer algebra system Mathematica. 
For Daubechies wavelets and symlets, = 6, . . . , 20, we find that Assump- 
tion 2.2 is indeed satisfied, and obtain accurate values of ct^ and f^, given 
in Table 1. 
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Algorithm 1 Verify assumption and compute constants 



/ ^ [0, 1] 
repeat 

calculate approximations 99^"^ to (^^"^ on /, n = 0, 1, 2 
deduce bounds on (c^)^"^ on /, n = 0, 1, 2 
deduce bounds on ct^ and v^p 

I smallest interval known to contain all maxima of cr^ 

i ^ i + 1 

until desired accuracy reached 

if a'^ bounded below zero on / then 

Assumption 2.2 is verified 
end if 



Daubechies Symlet 



N 




4 








4 




Vp 


6 


1, 


.251 


716 


0. 


,221 


993 


1 


.361 


961 


0, 


,106 


518 


7 


1, 


.276 


330 


0. 


,197 


328 


1 


.253 


835 


0, 


,248 


681 


8 


1, 


.250 


928 


0. 


,266 


316 


1 


.286 


722 


0, 


,173 


642 


9 


1, 


.222 


637 


0. 


,275 


519 


1 


.232 


334 


0, 


,302 


351 


10 


1, 


,199 


772 


0. 


,391 


629 


1 


.243 


114 


0, 


,255 


337 


11 


1, 


,195 


384 


0. 


,415 


019 


1 


.209 


007 


0, 


,324 


200 


12 


1, 


,189 


984 


0. 


,445 


388 


1 


.215 


480 


0, 


,335 


022 


13 


1, 


,182 


351 


0. 


,460 


792 


1 


.195 


567 


0, 


,385 


147 


14 


1, 


,172 


690 


0. 


,510 


179 


1 


.195 


969 


0, 


,405 


884 


15 


1, 


,165 


335 


0. 


,553 


767 


1 


.184 


307 


0, 


,446 


419 


16 


1, 


,159 


678 


0. 


,594 


027 


1 


.181 


901 


0, 


,465 


670 


17 


1, 


,154 


955 


0. 


,621 


941 


1 


.174 


105 


0, 


,496 


485 


18 


1, 


,150 


103 


0. 


,652 


913 


1 


.170 


871 


0, 


,520 


228 


19 


1, 


,145 


393 


0. 


,686 


434 


1 


.164 


974 


0, 


,551 


765 


20 


1, 


,141 


050 


0. 


,722 


113 


1 


.161 


837 


0, 


,571 


150 



Table 1: Computed values of constants 
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A Proofs 

We will need the following result, which is a version of Theorem 1 in Hiisler 
(1999). The result concerns the maxima of centred Gaussian processes whose 
variance functions are periodic; such processes are called cy do stationary. In 
Hiisler's original result, the maxima of a sequence of processes was shown 
to converge to a Gumbel random variable. In our result, we will specialise 
to a single process, and show this convergence occurs uniformly. 

Lemma A.l. Let T = T{n) — )■ cxd as n —)■ oo. In the notation of Hiisler 
(1999), let (A1)-(A3) and (B1)'(B4) hold, for a fixed process X„(t) = X{t), 
not depending on n. Further let a = (3, and let Hiisler's condition (1) hold. 
Define 

u{t) = (T„/z"^(r/mr). 
Then for any tq > 0, we have 

nMniT) > ujr)) ^ 
1 - e-^ 

as n ^ oo. 

Proof. Our argument proceeds as in the proof of Theorem 1 in Hiisler (1999). 
Without loss of generality, we may assume that an = 1- For r < tq, u{t) > 
u{to) —7- oo, and by definition 

mTlJ'{u{T)) = T. 

The approximation errors in parts (i) and (ii) of Hiisler's proof are thus 
0{g{S)T) and 0{pct) respectively. In part (iii), we note that 

uirf = 21og(r/r) - loglog(T/r) + 0(1), 

so Hiisler's term (4) is of order 

ri+''(r/r)i+''u(r)2/" exp{-n(r)V(l + t)} 

= ^1+^ exp{-[(l - 7)/(l + 7) - ??] log(r/r) + o(log(r/T))} = o(t), 

and term (5) is of order 

T2(r/T)2<5(r'')exp{-(2iog(r/T) - iogiog(r/T))(i - 5{t^)] 

= 0{TH{T^)\0g{T/T))=0{T). 



sup 

refO.rnl 
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In Hiisler's final display, we may thus write 

IP(M„(r) < u{t)) = exp{-(l + o(l))r} + o(t). 

As the process X{t) does not depend on n, the error in each of these 
approximations depends only on n = u{t), and the above limits hold as 
u — )• oo. (This can be seen from the precise form of the errors, as given 
in Piterbarg and Seleznjev, 1994, §3.1, and in Hiisler's proof.) Since u is 
decreasing in r, the limits are therefore uniform in r small. 

Consider the function 



f{x,y;T) := log 



1 — exp(— (1 + x)t) 



+ y 



defined on < r < tq, |x| < ^, \y\ < ^(1 — exp(— iro))/ro. The derivatives 

exp(-(l + x)r) 6f _ 1 
6y 



SI 
6x 



exp / Sy exp / 

are finite, and continuous in x, y and r, so by the mean value inequality, for 
n large, 



log 



/(O,0;r) + o(l 
log 



+ 0(1). 



As the above limits are uniform in r < tq, the result follows. 



□ 



We now apply this result to a cyclostationary process, composed of scal- 
ing functions ip, which we can use to model the variance of estimators f{jn)- 

Lemma A. 2. Define the cyclostationary Gaussian process 

I \ ^ /, J \ nr i-i-d. 



X{t) := a, 
For any 70 € (0,1), j. 



00, 



^{t-k)Zk, Zfc ~ iV(0,l). 



sup 

7e(0,7o] 
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-1ti 



sup \Xit)\>^ + b{jn) 
.te[0,2J"] «Unj 



as n 



00. 



Proof. For fixed 7, the result is a consequence of Theorem 2 in Gine and 
Nickl (2010); the statement uniform over (0, 70] follows, replacing Theorem 1 
of Hiisler (1999) in Gine and Nickl's proof with Lemma A.l. The conditions 
of Gine and Nickl's theorem are satisfied by Assumptions 2.1 and 2.2, as 
follows. 
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(i) X has almost-sure derivative 

^'W ■.= o:^^Y.^'{t-k)Zk, 

k€Z 

SO is continuous. X' is also the mean square derivative: 

h-^E[{X{t + h)- X{t) - hX'{t)f] 

= h-^ i^Pi^ -h-k)-^{t-k)- hif'it - k)f , 

kez 

which tends to as — )• 0, since the sum has finitely many non-zero 
terms. 

(ii) For i = 0, 1, define functions fi{x) := on [0,1], having wavelet 
expansions 

k j>J k 

in our wavelet basis on [0,1], for some J > jo, 2'^ > 6K. As ip is 
twice continuously differentiable, and if and have compact support, 
by Corollary 5.5.4 in Daubechies (1992), ■0 has at least two vanishing 
moments. Thus 

/3i,fc = (^\V',,fc) = 0, 

and 

k 

For t G [0,1], let v{t) denote the vector {(pj,k{t)) G 1^^^ so fi{t) = 
(a*,f(t)). Given s ^t,we have 

(a°,t;(5)) = l = (a°,r;(t)), 
{a\v{s)) = s^t = {a\v{t)), 

so the vectors v{s), v{t) are linearly independent. 
For s, t G M, define 

rx{s,t) := CoY[X{s),X{t)], al{t) := Var[X(t)] = rx{t,t). 

Then, if s,t G [-K,K], 

rx{s,t) = CT^^ ^ V(s - k)'p{t - k) 
fcez 
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so by Cauchy-Schwarz, 

rx{s,tf <alis)alit). 

If s,t G [k — K,k + K] for some A; G Z, the same applies by cyclo- 
stationarity. If not, then as ip is supported on [1 — K,K], we have 
rxis,t) = 0. However, for any t £ [0, 1], (a\f(^ + 2""^t)) = 1, so 

al{t) = a-'2~'\\v{t)f>0, 

and by cyclostationarity the same holds for all t E M. We thus again 
obtain 

rx{s,tf <al{s)aj,{t). 

(iii) We have 

so by Assumption 2.2, supjg[o,i] ^x(^) ~ ^^^^ maximum is at- 

tained at a unique to G [0, !)• If to S (0, 1), this satisfies the conditions 
of the theorem directly; if not we may proceed as in Proposition 9 of 
Gine and Nickl (2010). o"^ is twice differentiable, 

2a^c7'xito) = (a2)'(to)a2(to)"i/2 = 0, 

and 

= {alnto)al{to)-'/^ < 0. 
Finally, let v'{t) denote the vector G M^. Then for t G [0, 1], 

{a\v'it)) = m = l, 

so 

E[X'{tof] = a-^Y.^'{to-kr (A.l) 

kez 

= a-^2-''\\v'C^ + 2-'to)f>0. 

(iv) Since (/9 has support [1 — K,K], 

sup \rx{s, t)| = 0. □ 

s,t:\s-t\>2K-l 

We may now bound the variance of f{jn)- We will show that the variance 
process is distributed as the process X from the above lemma, so can be 
controlled similarly. 
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Proof of Theorem 2.3. Let I„ := [0,2-'"]. The process 

Xn{t) := /(J'") -{(J") (2-Jn^), t e In, 

C{jn) 

is distributed as 



for Zj^fc A^(0, 1), so by an orthogonal change of basis, as 



In case (i), X„ is distributed as the process X from Lemma A. 2, so we are 
done. 

In case (ii), set J„ := [2^,2-^" — 2K], and Kn '■= In \ Jn- On J„, X„ is 
distributed as the process X from Lemma A. 2, and we have 

sup|X„(t)| > n ) < P ( sup|X„(t)| > u 

< P ( sup|X„(t)| >u] + ¥{ sup|X„(t)| > n ) , 
\teJn J \teKn 



so for UnUn) ■= x{'yn)/a{jn) + b{jn), 



P ( SUp|X„(t)| > UnUn) - P SUp|X(t)| > Un{jn) 



<P sup|X„(t)| >n„(j„) 
< 8K{1 - $(Cn„(j„))) 

with a constant C > depending on (p. This term is o(7„), so the result 
follows by Lemma A. 2, applied to the process X on J„. □ 



B Source code 



The following program implements Algorithm 1 in Mathematica 8 or above. 
Note that we bound v^p by bounding the numerator and denominator of (2.3) 
separately over /. By (A.l), the numerator is positive; to bound inside 
(0, oo), we must therefore bound cr^ below zero. To verify Assumption 2.2, 
it is thus sufficient that we establish a finite positive value of v^. 
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Main: : usage = 

"Main[w, p, d, jmax, kmax] = bounds on the parameters sigma'"2_phi and\n" <> 
"upsllon_phi for the Wavelet w, performing computations using p digits\n" <> 
"of accuracy, stopping once accurate to d digits, after jmax steps, \n" <> 
"after evaluating phi at 2"kmax locations simultaneously, or after\n" <> 
"the results become indeterminate due to lack of precision. \n\n" <> 
"Main[DaubechiesWavelet[6] , 200, 6, 100, 12] = \n" <> 
" {Interval [{1.251716, 1.251716}], Interval [{0.221993, 0.221993}]}" 
Main [wavelet. , precision., digits., jmax_, kmax_] := 

Cascade [Filter [wavelet , precision], 2, PhiParams [digits, jmax, kmax]] 

Filter: : usage = 

"Filter [w] = the cascade filter of a Wavelet w." 
Filter [w_, p_ :MachinePrecision] := 

2 Transpose [WaveletFilterCoef ficients [w, WorkingPrecision -> p]] [[2]] 

PhiParams :: usage = 

"PhiParams [d, jm, km] = a callback for use with Cascade, whlch\n" <> 

"computes the parameters sigma"2_phi and upsilon_phi to d d.p., with\n" <> 

"resource limits jm, km." 
PhiParams [d_ , jm_, km_] := 

Params [Function [{w, j, offset, phi}, 

Module [{PhiMN , s2, s2max, ssp, sp2, sspp, s22b, s22a, ba, zoom}, 
PhiMN[m_, n_] := Plus m If[m == n, phi [ [n+1] ] "2, phi[[m+l]] phi[[n+l]]]; 

s2 = PhiMN [0, 0]; ssp = PhiMN [0, 1]; 

s2max = IntervalMax [s2] ; 

sp2 = PhiMN [1, 1]; sspp = PhiMN [0, 2]; 

s22b = IntervalRange [sp2] ; 

s22a = -IntervalRange [sspp+sp2] ; 

ba = s22b / s22a; 

zoom = Thread [(!TrueQ[# < s2max] fe /a s2) && (!TrueQ[# != 0] & /a ssp)]; 
{{s2max, ba}, zoom}]], d, jm, km] 

Params: : usage = 

"Params [Callback, digits, jmax, kmax] = a callback function for\n" <> 
"use with Cascade. The function will pass its arguments to Callback, \n" <> 
"which should return estimates of parameters, and an array of Indices k\n"<> 
"to restrict to. The function will terminate the cascade algorithm if\n" <> 
"the returned parameters are accurate to digits d.p., the resource\n"<> 
"limits jmax, kmax are reached, or the computation is indeterminate . \n" <> 
"Otherwise, it will restrict the algorithm to an interval I containing\n" <> 
"the requested indices , and continue . " 
Params [Callback. , digits., jmax., kmax.] := 

Function[{w, j, offset, phi}, Module[{ret, zoom, len, i, 1}, 
{ret, zoom} = Callback [w, j, offset, phi]; 
len = Dimensions [phi] [ [3] ] ; 
Which [ 

MemberQ[ret, Indeterminate, Infinity], 

Print ["indeterminate, j = ", j] ; 

ret = If [MemberQ [#, Indeterminate, Infinity], 

Interval [{-Infinity, Infinity}], #] & /@ ret; 

Prepend[ret, True] , 
(And m (IntervalAccurate[#, digits] &) /a ret) I I 
(j == jmax tt (Print ["hit jmax"]; True)) I I 
(len >= 2"kmax && (Print["hit kmax, j = ", j] ; True)), 

Prepend [ret , True] , 
len == 2"j , 

{1, i} = LongestSubsequence [Join [zoom, zoom], # == False &] ; 
If[l < 2-(j-l), {False}, Prepend [Mod [#-1, 2-j]+l & /a {i+1, i-1}. False]], 
True, 
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Prepend [Through [{Min , Max}[Position[zoom, True]]], False]]]] 

IntervalAccurate: : usage = 

" IntervalAccurate [i , d] = True if the Interval i is accurate to d d.p." 
IntervalAccurate [int_, digits.] := 

Equal as Round [Through [{Min, Max}[int]], 10' (-digits)] 

Int ervalRange : : us age = 

"IntervalRange [is] = an Interval giving the range of the Intervals is." 
IntervalRange [ints_] := 

Interval 8 Through [{Min, Max}[ints]] 

Int er valMax :: usage = 

"IntervalMax[is] = an Interval giving the maximum of the Intervals is." 
IntervalMax [ints_] := 

Interval S Map [Max, Transpose [ints /. Interval -> Identity]] 

LongestSubsequence: : usage = 

"LongestSubsequence [x, p] = {1, i}, where 1 is the length of the\n" <> 

"longest consecutive subsequence of x whose members satisfy p, and i\n" <> 

"is the index of its first member." 
LongestSubsequence [x_ , crit_] := 

Block [{$RecursionLimit = Infinity}, Module [{m, i, 1}, 

m = LengthWhile[x, crit] ; 

If [m == Length [x], {m, 1}, 

{1, i} = LongestSubsequence [Drop [x, m + 1] , crit]; 
If [1 > m, {1, i + m + 1}, {m, 1}]]]] 

Cascade: : usage = 

"Cascade[w, n. Callback] runs the cascade algorithm with filter w,\n" <> 
"boimding f and its first n derivatives .\n\n" <> 

"The function Callback [w, j, offset, phi] should return a list ret.\n" <> 

"If First[ret] is True, cascade will halt, returning Rest [ret] . \n" <> 
"If First [ret] is False, the algorithm will continue, and if\n" <> 
"Rest [ret] = {a, b} is given, future evaluation of f will be\n" <> 
"restricted to the indices a, . . . ,b of phi.\n" 
Cascade [w_, n_, Callback.] := 

Module [{of f set , u, j, k, g, f, fs, eps, phi, a, b, pts}, 
offset = 0; j = 0; k = Length[w]/2; 

u = BuildU[w, n+1] ; g = InitG[k, n+1] ; fs = {BuildF[g]}; 
While [True , 

g = StepG[g, u] ; f = BuildF[g]; eps = BoundError [f s , g, u] ; 
phi = ApplyError [f , eps, k] ; AppendTo[fs, f ] ; offset *= 2; j += 1; 
ret = Callback[w, j, offset, phi]; If [First [ret] , Return [Rest [ret] ]] ; 
{a, b> = If [Length [ret] > 1, Rest [ret], {1, Dimensions [phi] [ [3] ]}] ; 
If [a > b && Dimensions [phi] [ [3] ] == 2"j, 

{fs, g} = WrapFG[fs, g, k] ; b = b + 2* j ; offset -= 2"j] ; 
g = Take[g, All, All, {a, 2(k-l)+b}]; 
fs = RestrictF[fs, offset, a, b, k] ; offset += a-1]] 

BuildU: : usage = 

"BuildU[u-(0), n] = u-(0:n)" 
BuildU [w_, n_] := 

Map [Function [un, Table [un[ [m; ; All ; ; 2] ] , {m, 2}]], 

With[{signs = Table[(-l)"i, {i. Length [w] }]} , 

NestList[2 signs Accumulate [signs #] ft, w, n]]] 

InitG: : usage = 

"InitG[k, n] = g_0-(0:n)" 
InitG [k_, n_] := 

Constant Array [Reverse @ IdentityMatrix[2k-l] , n+1] 



15 



StepG: : usage = 

"StepG[g_j~(0:n), u-(0:n)] = g_{j+l}-(0:n)" 
StepG [g_, u_] := 

MapThread[Function[{gn, un}. Function [gni. 

Riffle S9 (ListConvolve[#, gni] ft) /a un] /S gn] , {g, u}] 

BuildF: : usage = 

"BuildF[g_j-(0:n)] = f_j-(0:n)" 
BuildF [g_] := 

Maplndexed [Function [{gn, part}. Module [{n, bn}, 

n = First [part] -1; bn = Table [Binomial [n, i](-l)"i, {i, 0, n}] ; 

Transpose [ListConvolve[bn, #, 1, 0] & /8 Transpose 8 gn]]], g] 

BoundAlpha: : usage = 

"BoundAlpha[g_j"(0:n), j] = alpha_j " (0 :n-l) " 
BoundAlpha [g_ , j _] : = 

l-Log[2, Max m Plus @S Abs S #] /j & /@ Rest S g 

BoundC: : usage = 

"BoundC[f_{0: j-l}-(0:n) , u-(0:n), alpha_j*(0:n-l)] = C_j-(0:n-l)" 
BoundC [f_, u_, alpha_] := 

MapThread [Function [{fn, un, an}, If[an <= 0, Infinity, 
(Max a Maplndexed[2-((an-l) (First[#2]-1)) Abs[#l] &, f n] ) 
(Max m Plus m (Abs [Accumulate [#]-!] &) /S un)/(l-2-(-an))]] , 
{Rest 9 Transpose 9 f. Most 9 u, alpha}] 

BoundError : : usage = 

"BoundError[f_{0:j-l}-(0:n), g_j*(0:n), u-(0:n)] = eps.j* (0:n-l) " 
BoundError [f_, g_, u_] := Module [{j, alpha, c}, 

j = Length [f]; alpha = BoundAlpha [g, j] ; c = BoundC [f, u, alpha]; 

MapThread[Function[{cn, an}, cn 2"(-j an)], {c, alpha}]] 

ApplyError: : usage = 

"ApplyError [f _j " (0 :n) , eps_j " (0 :n-l) , k] = intervals bounding f_j"(0:n-l)" 
ApplyError [f_ , eps_, k_] := 

MapThread [Function [{fn, en}, 

Map [Interval [{#-en, #+en}] &, fn, {2}]], 

{Most Take[f, All, All, {2k-l, -1}], eps}] 

WrapFG: : usage = 

"WrapFG[f_{0: j}"(0:n) , g_j*(0:n), k] = f and g wrapped around at integers" 
WrapFG [f_, g_, k_] := 

Module [{Wrap} , Wrap [x_] : = 

Map [Function [xn. Map [Function [xni, Join[xni [[1] ] , xni [[2] ][ [2k-l ;;]]]] , 
Partition [xn, 2, 1, {-1, 1}, {Constant Array [0 , Dimensions [x] [[3] ]]}]]] , 
x] ; {Wrap /9 f , Wrap 9 g}] 

RestrictF: : usage = 

"RestrictF[f_{0:j}-(0:n), offset, a, b, k] = f_{0: j}-(0:n) |offset+[a, b] " 
RestrictF [fs_, offset., a_, b_, k_] := 
Maplndexed [Function [{f j , part} , 

Module [{scale, start, end}, scale = 2" (Length [fs] -First [part] ) ; 
start = Floor [(of fset+a-l)/scale] -Floor [of f set/scale] +1 ; 
end = 2(k-l)+Floor[(off set+b-l)/scale]-Floor[offset/scale]+l; 
Take[fj, All, All, {start, end}]]], fs] 
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